Projeto final

Libraries

Stocks gathering

my_stocks <- c("ITUB3.SA", "ENBR3.SA",
               "FLRY3.SA", "B3SA3.SA", "BBDC3.SA",
               "EGIE3.SA", "TAEE11.SA", "WEGE3.SA", 
               "PSSA3.SA", "^BVSP")

all_asset_returns <- 
  BatchGetSymbols(my_stocks,
                first.date = "2000-01-01",
                bench.ticker = "^BVSP",
                type.return = "log",
                freq.data = "daily"
                )$df.tickers %>% 
  select(ref.date, ticker, ret.closing.prices, price.close) %>% 
  drop_na()

tickers_filter_table = tibble(
      stock = all_asset_returns$ticker %>% 
        unique,
      add_date = c(0, 2300, 2000, 0, 3000, 4000, 0)
    )


do_filter_ref_date <- function(stock, filter_table) {
  stock %>% 
      filter(ref.date > (
        stock$ref.date[1] + 
          filter_table$add_date[filter_table$stock == stock$ticker[1]]
        )
      )
}

asset_returns <- read_csv("\n", col_names = names(all_asset_returns)) %>% 
  mutate(
    ref.date = as.Date(ref.date),
    ret.closing.prices = as.numeric(),
    price.close = as.numeric(),
  )


for (t in unique(all_asset_returns$ticker)) {
  asset_returns <- asset_returns %>% 
    add_row(all_asset_returns %>% 
              filter(ticker == t) %>% 
              drop_na('ret.closing.prices') %>% 
              do_filter_ref_date(tickers_filter_table))
}

Funções - retornos e acf

plot_returns <- function(stock) {
  if (nrow(stock) > 0) {
    df_plot <- stock
    
    par(mfrow=c(1, 3))
    
    ts.plot(df_plot$price.close, main = paste("Ticker", stock$ticker[1], " closing prices"))
    ts.plot(
      df_plot$ret.closing.prices, 
      main = paste("Ticker", stock$ticker[1], "ret")
            )
    ts.plot(
      df_plot$ret.closing.prices^2,
      main = paste("Ticker", stock$ticker[1], "ret^2")
    )
  }
}

plot_asset_returns <- function (asset_returns) {
  tickers <- asset_returns$ticker %>%
    unique

  for (t in tickers) {
      asset_returns %>% 
        filter(ticker == t) %>% 
        plot_returns()
    }
}

box_test <- function (df, lag = 1, type = 'Ljung-Box') Box.test(df, lag = lag, type = type)

descriptive_asset_evaluation <- function (asset_returns) {
  tickers <- asset_returns$ticker %>% 
    unique
  
  for (t in tickers) {
    df_box <- asset_returns %>% 
      filter(ticker == t)
    
    acf(df_box$ret.closing.prices, main = paste("ACF", t)) %>% 
      autoplot() %>% 
      theme(
        plot.title = element_text(vjust = -50)
      )
    
    print(box_test(df_box$ret.closing.prices))
    print(adf.test(df_box$ret.closing.prices))
    print(ArchTest(df_box$ret.closing.prices))
  }
}


visualize_residual_distribution <- function(asset_returns) {
  tickers <- asset_returns$ticker %>% 
    unique
  
  par(mfrow = c(1,2))
  for (t in tickers) {
    ret <- (asset_returns %>% 
        filter(ticker == t) %>% 
        drop_na('ret.closing.prices'))$ret.closing.prices
    
    h <- hist(ret, 
              breaks=20, 
              col="red", 
              xlab="", 
              main= paste("Histogram", t)
          ) 
    xfit <- seq(min(ret),
                max(ret),
                length = 40) 
    yfit <- dnorm(xfit,
                  mean = mean(ret),
                  sd = sd(ret)) 
    yfit <- yfit * diff(h$mids[1:2]) * length(ret) 
    lines(xfit, yfit, col="blue", lwd=2)
    
    qqnorm(ret, pch = 1, frame = FALSE)
    qqline(ret, col = "steelblue", lwd = 2)
    
    print(paste("Realizando teste de normalidade para o ticker", t))
    print(shapiro.test(as.vector(ret[c(1:4999)])))
  }
  
}

calculate_assets_beta <- function (asset_returns) {
  asset_returns %>%
    group_by(ticker) %>%
    do(model = lm(ret.closing.prices ~ market_ret$ret.closing.prices,
              data = .))
}

Séries de fechamento, retornos e retornos quadrados

plot_asset_returns(asset_returns)

Os retornos aparentemente são estacionários. Iremos fazer o teste de ljung-box para confirmar o ruído branco nos resíduos e também realizar o teste de raíz unitária.

ACF

descriptive_asset_evaluation(asset_returns)

## 
##  Box-Ljung test
## 
## data:  df
## X-squared = 8.4382, df = 1, p-value = 0.003674
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_box$ret.closing.prices
## Dickey-Fuller = -17.489, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  df_box$ret.closing.prices
## Chi-squared = 273.49, df = 12, p-value < 2.2e-16

## 
##  Box-Ljung test
## 
## data:  df
## X-squared = 40.974, df = 1, p-value = 1.543e-10
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_box$ret.closing.prices
## Dickey-Fuller = -16.892, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  df_box$ret.closing.prices
## Chi-squared = 527.07, df = 12, p-value < 2.2e-16

## 
##  Box-Ljung test
## 
## data:  df
## X-squared = 3.0717, df = 1, p-value = 0.07967
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_box$ret.closing.prices
## Dickey-Fuller = -14.665, Lag order = 16, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  df_box$ret.closing.prices
## Chi-squared = 824.95, df = 12, p-value < 2.2e-16

## 
##  Box-Ljung test
## 
## data:  df
## X-squared = 23.779, df = 1, p-value = 1.081e-06
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_box$ret.closing.prices
## Dickey-Fuller = -17.779, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  df_box$ret.closing.prices
## Chi-squared = 387.87, df = 12, p-value < 2.2e-16

## 
##  Box-Ljung test
## 
## data:  df
## X-squared = 3.2264, df = 1, p-value = 0.07246
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_box$ret.closing.prices
## Dickey-Fuller = -14.784, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  df_box$ret.closing.prices
## Chi-squared = 642.15, df = 12, p-value < 2.2e-16

## 
##  Box-Ljung test
## 
## data:  df
## X-squared = 2.8955, df = 1, p-value = 0.08883
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_box$ret.closing.prices
## Dickey-Fuller = -10.607, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  df_box$ret.closing.prices
## Chi-squared = 128.7, df = 12, p-value < 2.2e-16

## 
##  Box-Ljung test
## 
## data:  df
## X-squared = 3.4689, df = 1, p-value = 0.06253
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df_box$ret.closing.prices
## Dickey-Fuller = -16.147, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  df_box$ret.closing.prices
## Chi-squared = 1531.6, df = 12, p-value < 2.2e-16

Residual distribution

visualize_residual_distribution(asset_returns)

Busca por melhor modelo a partir do tune

calculate_portfolio_covariability <- function(tib_returns) {
    cov(tib_returns$returns, tib_returns$ret.closing.prices) / var(tib_returns$ret.closing.prices)
}

garchmodels que nao funciona

fit_ts_stock <- function (stock) {
    if (nrow(stock) > 0) {
      stock <- stock %>%
        timetk::future_frame(
          .length_out = 3,
          .date_var = ref.date,
          .bind_data = TRUE
        ) %>%
        mutate(ref.date = as.POSIXct(ref.date))
  
      stock_train <- stock %>%
          drop_na()
  
      stock_future <- stock %>%
        filter(is.na(ret.closing.prices))
  
  
      stock_model <- garchmodels::garch_reg(
        mode = "regression",
        arch_order = tune::tune(),
        garch_order = tune::tune(),
        ma_order = tune::tune(),
        ar_order = tune::tune(),
        tune_by = "sigmaFor"
      ) %>%
        parsnip::set_engine("rugarch")
  
      stock_recipe <- recipes::recipe(
        ret.closing.prices ~ ref.date,
        data = stock_train
      )
  
      stock_wflw <- workflow() %>%
        add_recipe(stock_recipe) %>%
        add_model(stock_model)
  
      stock_resamples <-  time_series_cv(
          stock_train,
          date_var = ref.date,
          initial = "1 year",
          assess = "3 months",
          skip = "1 month",
          cumulative = TRUE,
          slice_limit = 5
      )
  
      stock_tune_results <- tune_grid(
        object = stock_wflw,
        resamples = stock_resamples,
        param_info = dials::parameters(stock_wflw),
        grid = 5,
        control = control_grid(
          verbose = TRUE,
          allow_par = TRUE,
          parallel_over = "everything"
        )
      )
  
  
      stock_tune_results %>%
        tune::show_best(metric = "rmse")
    }
} 



results <- asset_returns %>% 
  group_split(ticker) %>% 
  map(fit_ts_stock)


results

Analisando CAPM

Qual o risco de cada ação em relação ao mercado?

evaluate_assets <- function(tickers, 
                            bench_ticker = "^BVSP", 
                            market = "^BVSP",
                            first.date = "2000-01-01", 
                            last.date = Sys.Date(),
                            assets_weights = NULL
                            ) {
  asset_returns <- BatchGetSymbols(tickers,
                first.date = first.date,
                last.date = last.date,
                bench.ticker = bench_ticker,
                type.return = "log",
                freq.data = "daily"
                )$df.tickers %>% 
  select(ref.date, ticker, ret.closing.prices) %>% 
    drop_na()
  
  tickers <- asset_returns$ticker %>% 
    unique
    
  for (t in tickers) {
    asset_returns %>%
      filter(ticker == t) %>%
      select(ref.date, ret.closing.prices) %>%
      drop_na()
  }
  
  
  portfolio_ret <- asset_returns %>%
  tq_portfolio(assets_col  = ticker,
               returns_col = ret.closing.prices,
               weights     = assets_weights,
               col_rename  = "returns")
  
  market_ret <-
    BatchGetSymbols(market,
               first.date = first.date,
               last.date = last.date,
               freq.data = "daily",
               type.return = "log")$df.tickers %>%
    select(ref.date, ret.closing.prices) %>% 
    as_tibble()
  
  ts.plot(portfolio_ret$returns)
  ts.plot(market_ret$ret.closing.prices)
  
  portfolio_market_ts <- portfolio_ret %>% 
     inner_join(market_ret, by = 'ref.date') 

  calculate_portfolio_covariability(portfolio_market_ts)
  
  betas <- asset_returns %>%
    group_by(ticker) %>%
    inner_join(market_ret, by = 'ref.date') %>% 
    do(model = lm(ret.closing.prices.x ~ ret.closing.prices.y,
              data = .))

  betas
}


assets_lm <- evaluate_assets(my_stocks)

assets_lm$betas <- map(assets_lm$model, coef) %>% 
  map_dbl(2)

assets_lm

## para calcular o beta do portfolio precisamos dos weights de cada uma das ações, e então multiplicar cada beta com o seu peso respectivo - leia-se beta_porfolio_byhand do arquivo capm_model.rmd

Dynamic model - CAPM

dlm_routine <- function (df, ticker_string, index_string) {
  my_dlm = function(parm, x.mat) {
    parm = exp(parm)
    return (
      dlmModReg(
        X = x.mat, 
        dV = parm[1], 
        dW = c(parm[2], parm[3])
      )
    )
  }
  

  ticker_ret <- df %>%
    dplyr::select('ref.date','ticker','ret.closing.prices') %>%
    dplyr::filter(ticker == ticker_string) %>%
    drop_na('ret.closing.prices')
  
  index_ret <- df %>%
    dplyr::select('ref.date', 'ticker','ret.closing.prices') %>%
    dplyr::filter(ticker == index_string) %>%
    drop_na('ret.closing.prices')
  

  rang = range(
    index_ret$ret.closing.prices,
    ticker_ret$ret.closing.prices
  )
  
  tib_returns <- ticker_ret %>% 
    inner_join(index_ret, by = 'ref.date')
  
  
  print(tib_returns)

  plot(
    tib_returns$ret.closing.prices.y,
    tib_returns$ret.closing.prices.x,
    xlab = paste0("Market return", index_string),
    ylab=ticker_string,
    xlim=rang,
    ylim=rang
  )
  
  capm_ticker <- lm(tib_returns$ret.closing.prices.x ~ tib_returns$ret.closing.prices.y)
  
  abline(
     capm_ticker$coef,
     col = 2,
     lwd = 3
   )
  
 title(
   paste(
     ticker_string,
     "=",
     round(capm_ticker$coef[1],4)," + ",
     round(capm_ticker$coef[2],4),
     "*SP500",
     sep=""
   )
 )
  
    
  fit <- dlmMLE(
    y = tib_returns$ret.closing.prices.x,
    parm = c(1, 1, 1),
    x.mat = tib_returns$ret.closing.prices.y,
    build = my_dlm,
    hessian = T
  )

  se = sqrt(exp(fit$par))

  mod_std  = my_dlm(fit$par, tib_returns$ret.closing.prices.y)
  mod_filt = dlmFilter(tib_returns$ret.closing.prices.x, mod_std)
  mod_smot = dlmSmooth(mod_filt)
  
  # colocamos o codigo original (com date -1) e não funcionou
  # colocamos sem e parece estar tudo ok 
  date = tib_returns$ref.date
  
  print(paste("isso", length(date)))
  print(paste("esse", length(mod_filt$m[,1][-1])))
  plot(
    date,
    mod_filt$m[,1][-1],
    xlab = "day",
    ylab = expression(alpha[t]),
    type = "l",
    main = ""
  )

  lines(
    date,
    mod_smot$s[,1][-1],
    col = 2
  )
  abline(
    h = capm_ticker$coef[1],
    col=3
   )
  abline(h = 1, lty = 2)
  
  
  plot(date,
       mod_filt$m[,2][-1],
       xlab = "day",
       ylab = expression(beta[t]),
       type = "l",
       main = ""
 )
  lines(date,
        mod_smot$s[,2][-1],
        col=2
  )
  
  abline(
    h=capm_ticker$coef[2],
    col=3
  )
  abline(h = 1, lty = 2)
}

dlm_routine(asset_returns, ticker_string = "ITUB3.SA", index_string = "^BVSP")

Evaluating garch

evaluate_garch <- function(ticker, data, mean.model, variance.model) {
  print(paste("Ticker", ticker, "usando os modelos"))
  print(mean.model)
  print(variance.model)
  garch_spec <- ugarchspec(
    mean.model = mean.model,
    variance.model = variance.model,
    distribution.model = "std"
  ) 
  
  garch_fit <- ugarchfit(spec = garch_spec, data = data$returns[-1])
  
  print(paste("Visualizando erros do modelo para o ticker", ticker))
  print(garch_fit@fit$coef)
  tibble(ticker = ticker, fit = c(garch_fit), AIC=list(as.vector(infocriteria(garch_fit))))
}

garch_forecast <- function(garch_fit) {
  garch_forecast <- ugarchforecast(garch_fit, n.ahead = 1)
  garch_forecast %>% 
    plot
}

make_asset_garchs <- function(asset_returns, garchOrders) {
  if (nrow(garchOrders) > 0 & nrow(asset_returns) > 0) {
    tickers <- asset_returns$ticker %>% 
      unique
    
    models <- tibble(ticker = NA, fit = NA, AIC = NA)
    
    for (t in tickers) {
      print(paste("Rodando garch para a ação ", t))
      ret <- (asset_returns %>% 
        filter(ticker == t) %>% 
          drop_na('ret.closing.prices') %>% 
          mutate(
            returns = ret.closing.prices
          ))
      
      mean.model <- (garchOrders %>% 
        filter(ticker == t))$mean.model
      
      variance.model <- (garchOrders %>% 
        filter(ticker == t))$variance.model
      
      fit_result <- evaluate_garch(t, ret, mean.model, variance.model)
      
      models <- models %>% 
        add_row(fit_result)
    }
    
    models
  } else {
    print("Assets e ordens dos garchs são necessários")
  }
}


garchOrders <- tibble(
  ticker = asset_returns$ticker %>% 
    unique(),
  mean.model = list(
    c(0, 0),
    c(0, 0),
    c(0, 0),
    c(0, 0),
    c(0, 0),
    c(0, 0),
    c(0, 0)
  ),
  variance.model = list(
    c(1, 2),
    c(1, 2),
    c(1, 2),
    c(1, 2),
    c(1, 2),
    c(1, 2),
    c(1, 2)
  )
)

result <- make_asset_garchs(asset_returns, garchOrders)

for (t in unique(asset_returns$ticker)) {
  print(paste("Info criteria do ticker", t))
  print(infocriteria((result %>% 
    filter(ticker == t))$fit[[1]]))
}